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It is well known from the quantum theory of strongly correlated systems that poles (or more 
subtle singularities) of dynamic correlation functions in complex plane usually correspond to the 
collective or localized modes. Here we address singularities of velocity autocorrelation function Z 
in complex explain for the one-component particle system with isotropic pair potential. We have 
found that naive few poles picture fails to describe analytical structure of Z(cj) of Lennard-Jones 
particle system in complex plain. Instead of few isolated poles we see the singularity manifold 
of Z(uS) forming branch cuts that suggests Lennard-Jones velocity autocorrelation function is a 
multiple-valued function of complex frequency. The brunch cuts are separated from the real axis 
by the well-defined “gap”. The gap edges extend approximately parallel to the real frequency axis. 

The singularity structure is very stable under increase of the temperature; we have found its trace 
at temperatures even several orders of magnitude higher than the melting point. Our working 
hypothesis that the branch cut origin is related to the “interference” in Z of one-particle kinetics 
and collective hydrodynamic motion. 

PACS numbers: 61.20.Ne, 65.20.De, 36.40.Qv 


I. INTRODUCTION 

Dynamic correlation functions (DCF) are one of the 
main tools that allow understanding nature of condensed 
matter particle systems W- Fourier spectra of DCF 
keep the information about the spectrum and inverse life¬ 
time of collective excitations, particle diffusion and most 
other key properties of the system. As a rule, using the 
results of numerical simulations, like molecular dynamics, 
one can find spectrum of DCF on real (or sometimes on 
imaginary) axis in the frequency cj-space. However, in¬ 
teresting fitches of the DCF should be hidden at complex 
uj. It is well known from the quantum theory of strongly 
correlated systems that poles (or more subtle singulari¬ 
ties) of DCF in complex plane usually correspond to the 
collective or localized modes. Then, for example, the 
real part of the pole position in the cj-plane produces the 
energy of the excitation while the imaginary part corre¬ 
sponds to the inverse life time mm- Interesting question 
what singularities of DCF(cj) for classical particle system 
one can find in the complex cj-plain. 

Here we consider the onecomponent particle system 
with isotropic Lennard-Jones (LJ) pair potential and 
focus mainly on the velocity autocorrelation function 
(VAF) Z{t). It is well known that in liquid phase Z(t) is 
nonmonotonic at short time scales t ~ tq, where ro is the 
period of the particle motion in the effective potential 
well formed by the surrounding particles (i.e. the in¬ 
verse Einstein frequency) [2 J 0 6 . In the dilute gas phase 
and in the supercritical fluid far above melting and criti¬ 
cal temperatures Z(t) decays monotonically with time at 
time scale of the order of the relaxation time of the par¬ 
ticle diffusion: r = mD / k B T , where D is the diffusion 
coefficient, T is the temperature and m is the particle 


mass [2,|5j|6]. 

There are many approximations of Z(t) for simple par¬ 
ticle systems. It is well established that for satisfactory 
approximation of Z(t) one should take more than one 
relaxation time in the memory function or even the con¬ 
tinuum of the relaxation times mmm. But how then 
the manifold of relaxation times and Einstein frequencies 
look like? Here we search for answers to these questions. 

As far as we know there is no universal explicit ex¬ 
pression that produces Z(t) equally well at small and 
hydrodynamic (large) time scales [2j [6]. On the other 
hand, low approximation accuracy do not allow reliable 
investigation of singularity manifolds in the complex em¬ 
plane. Therefore here we investigate Z(uj) numerically 
and develop the machinery for numerical analytical ap¬ 
proximation. 

We see the singularity manifold of Z(co) forming branch 
cuts that suggests LJ VAF is a multiple-valued func¬ 
tion of complex frequency. The brunch cuts are approx¬ 
imately parallel to the real frequency axis and separated 
from it by the well-defined “gap”. The singularity man¬ 
ifold is stretched along the real axis. Going higher and 
higher with temperature the singularities more and more 
group and Z(uj) better and better agrees with the expo¬ 
nential memory function approximation of Z(t) [2- The 
singularity structure is very stable under increase of the 
temperature; we have found its trace at temperatures 
even several orders of magnitude higher than the melt¬ 
ing point. 

When t > th, where ty is some characteristic tran¬ 
sient time for hydrodynamic regime, Z(t) has nontrivial 
power law decaying tail oc t -3 / 2 . In frequency space that 
causes nonanalyticity of Z{uj) ~ cj -1 / 2 at small uj [TT] . 
The short-time t <ty behaviour of Z(t) with satisfactory 
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accuracy can be simulated by a number of damped and 
over damped oscillators that formally produce exponen¬ 
tial long-time time decay. In frequency representations 
these oscillators produce analytical function with a num¬ 
ber of poles. Binding “analytical” oscillators with “non- 
analytical” hydrodynamics at t ~ th produces, from our 
point of view, nonanalyticity in the Fourier transform of 
Z(t). We do see that the time scale th well corresponds 
the half-width Ao; of the gap between branch cuts in 
Z(uj) in wide temperature range: th ~ 2tt/Auj. 

II. MODEL CALCULATIONS 

One way to go into the complex cj-plain is the z- 
transform of DCF(t) [this is fast complex-u; Fourier trans¬ 
form]. This approach has been used in Refs. (TOj fl2] . 
However the analytical continuation of DCF were not 
there the purpose of the study except the answer to 
the question if the singularities of the memory func¬ 
tion belong to the stability manifold |z| < 1 of the z- 
transform. Since the stability of ^-transform is limited 
in the complex cj-plain one should search for alternatives. 
Other traditional methods of the analytical continuation, 
like integration of Cauchy-Riemann equations or differ¬ 
ent methods of series reexpansion m, also are unstable 
approaching DCF(cj) singularities. 

The promising way to study the singularities of DCF 
in the complex cj-plain is to built at real uj an approx¬ 
imation of DCF by a meromorphic function and finally 
do the analytical continuation of it. [In complex analy¬ 
sis, meromorphic function is a function that is holomor- 
phic except a set of isolated points [141 [T5] .] The Pade- 
approximation is the keystone of one of the analytical 
continuation methods that follows this receipt [T6MT9] . 
Here we perform complex-o; spectroscopic investigation 
of Z{uS) based on the Pade-approximation, while Z(t) we 
find from Molecular Dynamic (MD) simulations. 

For MD simulations of Z(t), we have used DL.POLY 
Molecular Simulation Package ED] developed at Dares- 
bury Laboratory. For simulations we use the LJ pair 
potential model in a wide range of parameters. For 
LJ liquid we apply the standard pair potential, U(r) = 
4£[(cr/r) 12 — (cr/r) 6 ], where 6 is the unit of energy, and 
a is the core diameter. In the remainder of this paper 
we use the dimensionless quantities: f == r/cr, U = U/e, 
temperature T = T/e, density p = 7Vcr 3 /U, and time 
t = t/[(Ty/m/e\, where m and V are the molecular mass 
and system volume correspondingly. As we will only use 
these reduced variables, we omit the tildes. 

For simulations, we have considered the system of 
N = 128 3 ~ 2.1 • 10 6 particles that were simulated un¬ 
der periodic boundary conditions in 3-dimensional cube 
mostly in the Nose-Hover (NVT) and also in NPT en- 
sambles. Such a large number of particles is necessary 
to correctly describe long time behaviour of correlation 
functions, see Ref. [21]. We consider the system at fixed 
density p — 1 and different temperatures in the interval 


T G (1.4,200). According to equilibrium temperature- 
density phase diagram [221426] , this range covers thermo¬ 
dynamic states from the liquid just above the melting 
line up to the supercritical fluid approaching the ideal 
gas limit. The MD time step was t = 0.0001 — 0.001 cho¬ 
sen so that to provide good energy conservation for given 
thermodynamic conditions. 

To obtain perfect hydrodynamic tails of Z(t) we im¬ 
proved the code of DL_POLY Molecular Simulation Pack¬ 
age [20] and inserted inside specially designed parallel 
MPI-code to calculate Z(t) for very large systems [2Tj . 
Calculations of Z{t) tails for 2 • 10 6 particles requires at 
least 128 processors with 10 —15 Gb of operational mem¬ 
ory per each one. 


III. RESULTS 

A. Fluid just above the melting line 

1. The branch cut 

We start our investigation from T = 1.4 and p = 1: 
these parameters correspond to fluid, just above the melt¬ 
ing line. The results are shown in Figs. 0-1 m 
obtained using MD simulation is shown in Fig. m a) and 
Fi g .0;b) represent Z{uj) obtained by Fourier transforma¬ 
tion of Z(t). At t > 1, Z(t) becomes positive and as it 
should be in fluid m, and at t > 5 it demonstrates t 3 / 2 
asymptotic behaviour, see insert in Fig.[lj[b). We use this 
long time asymptotic to obtain good Fourier transform 
of Z(t). We perform extrapolation of Z(t) to long times 
by at -3 / 2 asymptotic with the appropriate value of a. 
[Of course, we have tested that this procedure does not 
change Z(cj) behavior and does not influence the prop¬ 
erties of analytical continuation to complex cj-plain.] As 
the result we obtain smooth Z(uj) curve at all interesting 
values of uj (see Fig.JTJb)). The Z(uS) curve has the form 
typical to that for simple liquids 0©. In particular, it 
demonstrates pronounced maximum at uj = o; max . 

Figs. [T] (c) and (d) show analytical continuation of 
Z(yj) into complex cj-plain using multipoint Pade ap¬ 
proximation built on top of 1400 uniformly distributed 
knot-points in (0, 6.3o; max ). Technical aspects of Pade 
approximation can be found in Sec.[V] Building the Pade 
approximant we explicitly take into account that Z(uj) is 
even function at real uj. So the continued fraction of Pade 
approximant is the function of uj 2 . 

3D plot of \Z(uj)\ in the complex u;-plane is shown in 
Figs, [l] (c) and (d). Regular behaviour of Z(uj) for small 
imaginary frequencies ends abruptly by the “walls” of 
singularities constructed from poles (and zero nodes) of 
the Pade-approximant. Fig. IT] (d) shows plot of In \Z(uj)\ 
in the domain Rea;, Ima; G (—4.5a; max , 4.5a; max ). Such 
series of poles and zeros is the way the Pade approxi¬ 
mation typically represents brunch cuts of multi-valued 
functions (see Sec. |v|). 
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FIG. 1. VAF at real frequency axis (a) and in time representation (b) for T — 1.4 and p — 1: results of MD simulation. The inset 
in (b) shows long time tail of Z(t) in the double logarithmic scale; the red bullet • points the time scale At = 2 tt/Auj, where Acj is 
the half-width of the branch cut. Graphs (c) and (d) show analytical continuation of Z(uS) into complex cj-plain using multipoint 
Pade approximation built on top of 1400 uniformly distributed knot-points in (0, 6.3cj ma x)- 3D plot of \Z(uo)\ in the complex 
cj-plane; the black curve in (a) is Z(uS) at real cj. Regular behaviour of Z(uS) for small imaginary frequencies ends abruptly by 
the “walls” of singularities. Graph (d) shows In \ Z(uj)\ in the complex frequency domain Reca, Imo; G (—4.5 cj max 5 4.5cJ ma x) • 


In the insert in Fig. [TJb) we show the long-time be¬ 
haviour of Z(t) in double logarithmic scales. The red 
bullet points the time scale At = 27 r/Aco , where, we re¬ 
mind, Au: is the characteristic scale approximately equal 
to the half-width of the gap between branch cuts. As 
follows, At ~ th ,where ty corresponds to crossover of 
system dynamics from kinetic to hydrodynamic regime. 


2. Hydrodynamic asymptotic of VAF: additional branch cut 
at small frequencies 

As we have already mentioned, at timescales much 
larger than inverse Einstein frequency VAF is positive 
and it has the following asymptotic behaviour: Z{t) ~ 
t- 3 / 2 [2ED- It generates singular terms in Z(uj —» 0) ^ 
const — yj\io\ at real frequency axis, see Fig. pk for il¬ 
lustration. This nonanalyticity at uj = 0 should"" produce 




FIG. 2 . (a) Hydrodynamic tales of Z(t) ~ t 3 ^ 2 generate 

additional singular terms in Z(pj 0) ~ const — y/\oJ\ at 
real frequency axis. This nonanalyticity at uj — 0 produces 
the branch cut shown in (b) where we use the same Pade 
approximant as in Fig. [I] 


additional branch cut. We do see it in Fig. [2^ where Z(oj) 
is shown. Preparing Fig. It we have used the same Pade 
approximant as we have used working on Fig. Bc)-(d). 
Hydrodynamic branch cut is on the imaginary axis and 
so it lies transversely to the branch cut presented in the 
Figs. [ljc)(d). Note that for liquid near the melting line 
hydrodynamic singularity is located at only small vicin¬ 
ity of uj = 0 (compare Fig. 0a) and Fig. [2|a)). Thus 
the hydrodynamic branch cut is only detectable at small 
frequency scales (compare Fig. 0 c)(d) and Fig.0d)). 


B. From fluid to gas: evolution of Z(uf) 

Below we test how stable is the the branch cut in Z(uS) 
when we increase the temperature of the fluid far above 
the melting line. It follows that the branch cut is very 
stable to temperature. 


1. T — 40 and p — 1 

For density p — 1 the temperature T = 40 is charac¬ 
teristic temperature when fluid local structure vanishes, 
see Ref. [26]. However the branch cut in Z (uj) is still well 
observable, see Fig. [3j We see from Fig. |3j]c) that the 
“small” cut originating from a/cJ singularity at uj 0 
continuously transforms into the “large” branch cut that 
goes parallel to the real cj-axis. 

Insert in Fig. [3jb) shows Z(t) in double logarithmic 
scale. Dash-dotted line there sketches the hydrodynamic 
oc l/t 3 / 2 -tail. The red bullet shows At. It follows that 
again At th • 
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FIG. 3. VAF at real frequency axis (a) and in time rep¬ 
resentation (b) for T — 40 and p — 1: results of MD simu¬ 
lation. Insert in (b) shows Z(t) in double logarithmic scale. 
Graphs (c) and (d) show \Z(u)\ and ln|Z(cj)| analytically 
continued into complex c^-plain using multipoint Pade ap¬ 
proximation built on top of 1400 uniformly distributed knot- 
points in u G (0, 200). The (half) width of the branch cut Ac o 
we mark by the blue bullet • on Z(cj). Characteristic time 
At = 2n/Aio we show by the red bullet • on Z(t) curve. 


FIG. 4. VAF at real frequency axis (a) and in time rep¬ 
resentation (b) for T = 200 and p — 1: results of MD sim¬ 
ulation. Graphs (c) and (d) show 3D and density plots of 
\Z(u)\ in complex c^-plain. Multipoint Pade approximation 
has been built on top of 1400 uniformly distributed knot- 
points in uj G (0,500). Again, the (half) width of the branch 
cut Auj we mark by the blue bullet • on Z (co>). Characteristic 
time At = 2tt/Aco we show by the red bullet • on Z(t) curve. 
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2. T — 200 and p — 1 

For density p = 1 and temperature T = 200 we have 
the slightly nonideal gas. However even at such high 
temperatures there is a branch cut, very small one, as 
follows from Fig. |4j Insert in Fig.[4^b) shows Z(t) in dou¬ 
ble logarithmic scale, where the dash-dotted line sketches 
oc l/t 3 / 2 -tail. The red bullet we put at t = At. Again, 
At ~ t h . 

Except small branch cuts, Fig. |4jd,e) reveal two iso¬ 
lated poles located on the imaginary axis. The appear¬ 
ance of such poles at hight temperatures shows that the 
dynamics of the system is near to that for the ideal gas. 
Indeed the simplest low-density-limit exponential rela¬ 
tion for Z(t) obtaining from either the Enskog approx¬ 
imation or the Brownian one has the same analytical 
structure with two pure imaginary poles [2]. 


IV. DISCUSSION 


Velocity autocorrelation function in general can be ex¬ 
pressed as follows in the Fourier space [2]: 


Zip) 


A 

—iui + M(u>) ’ 


( 1 ) 


where A is constant and M{ui) is the “memory” function. 
Here “tilde” above Z means that we put Z(t < 0) = 0. 
If we say that Z depends on \t\ then Z(uS) = (Z(cj) + 

Z(~u))/2. 

Using the projector operator formalism ca© it is pos¬ 
sible to find an exact representation of the memory func¬ 
tion as the continued fraction of the form 


M(u) = 


Mi( 0) 


—iuj + 


m 2 (o) 5 

— ... 


( 2 ) 


where M n (cj), n > 0 is the hierarchy of memory func¬ 
tions. The coefficients M n ( 0) are related to the frequency 
moments and may be in principle calculated through in¬ 
teraction potential and static properties. In practice the 
only few first moments can be calculated and so one usu¬ 
ally has to truncate the continued fraction Q at some 
finite term [28]. The simplest case of the first-order trun¬ 
cation M n (0) = 0, n > 1 gives trivial exponential decay 
of VAF; the second one corresponds to non-trivial case 
of M{uj) oc 1 /(—iu) + 1 /t) [2] which demonstrates qual¬ 
itatively correct VAF behaviour but quantitatively fails 
to approximates Z(uS) even at real uj. As a matter of 
fact, no finite truncation scheme gives quantitative de¬ 
scription at whole co range. This leads researchers to use 
phenomenological approximations for memory functions, 
see for the review. The general conclusions about 
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such approximations is that one relaxation time models 
are not enough to describe Z(uj) satisfactory; at last two 
relaxation times is needed. Below we consider one of such 
approximation as well as alternative approach based on 
mode coupling theory and investigate what behaviour of 
Z in the complex uo plain these approaches produce. 

It should be noted here that the time evolution of Z(t) 
can be qualitatively understood if we truncate the con¬ 
tinued fraction of M. Then Z{uj) has few poles. When 
the poles of Z(u) are purely imaginary then Z(t) should 
decay monotonically and the poles correspond to the re¬ 
laxation times; the nonzero real part of the poles induce 
nonmonotonic behavior of Z(t) [2]. Unfortunately this 
approximation usually has very poor accuracy [2.. These 
considerations also fail explaining the hydrodynamic time 
scales, t r, where Z(t) shows the universal long-time 
tails Z(t) governed by hydrodynamic fluctuations, see 
Refs. mEUESMI]. 


A. VAF: Two-exponential approximation of the 
memory function 

There is well known approximation involving two re¬ 
laxation times, where [2 El]: 

M(t > 0) ~ At 4 e~ at + Be~ bt \ (3) 



00 1 0 20 Re(co)/co max 50 



FIG. 5. (a) The blue dotted curve is VAF obtained by MD 

simulation for T = 1.4 and p — 1 while the red curve is the ap¬ 
proximation when the memory function contains two distinct 
relaxation processes. Graphs (b) and (c) show the absolute 
value and the complex argument of the approximation ana¬ 
lytically continued to the complex plain. 


Taking A, B, a and b as an adjusting parameters and do¬ 
ing Fourier transform of M(t) we can fit VAF. The result 
is illustrated in Fig.[5ja) for T = 1.4 and p = 1. The fit is 
not very good however it seems from the first glance that 
main features of VAF this approximation reproduces, at 
least at moderate and large frequencies. However going 
to the complex cj-plain we see absolutely different be¬ 
haviour than exact VAF shows, see Fig. UN and (c): 
there is no branch cut parallel to the real axis. Below we 
investigate better approximation, however it also does 
not show coincidence with the exact result in the com¬ 
plex plain. 


B. Mode coupling approach and viscoelastic 
approximation 


An alternative way to calculate VAF based on mode 
coupling theory says that Z(t) can be expressed in the 
form : 

o roo 

Z{t)= 2^p J 0 f( k ) F s(.k,t)[Ci(k,t) + 2C t (k,t)]dk, 

(4) 


where F s is the self part of intermediate scattering func¬ 
tion, while Ci and Ct are the longitudinal and transverse 
current autocorrelation functions [2]. The wight-function 
f(k) is provides the appropriate cut-off of the integral at 
large enough k [33] . 

The analytical calculation of the integrand in @ is 
the non-trivial task. The only F s (k,t ) can be estimated 
relatively easy within the framework of the gaussian ap¬ 
proximation: F s (k,t ) ~ exp(— a(t)k 2 ), where a(t) is an 
unknown function which may be related to either the 
mean-square displacement [2] or VAF itself [7]- The cal¬ 
culation of Ci, C t is a more complicated task. Here we 
use the simplest approximation - viscoelastic one. 

For Ci we have used the following standard expres¬ 
sion [2j: 


1 uj 2 ujq Re Ni 

n (-W Im Ni — ^|y + 2 + (u Re AT ,) 2 ’ 


N t (u,k) 


-iuj + 7;(fc) 


(5) 

(6) 


Here coq = k \JTjm and S(k) is the structure factor that 
we take from MD simulation. The effective frequencies 
are defined as follows [ 33] : 


uf{k) = 3 ul(k)+ 

2 f 3sin (kr 0 ) | 6sin(fcr 0 ) 6cos(fcro) \ , . 

WE V kr 0 (fcr 0 ) 3 (fcr 0 ) 2 J ’ 

where uoe and tq are the adjusting parameters of the 
order of the Einstein frequency and interparticle spacing. 
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FIG. 6. The black curve in (a) is VAF obtained within MD simulation while the red curve in (a) shows VAF obtained 
within the viscoelastic model. Analytical continuation of VAF in viscoelastic model into complex frequencies: (b) represents 
the amplitude of VAF while (c) is its complex phase. No brunch cut is seen. Graphs (d) and (e) show behaviour of Z(cj) at 
small (hydrodynamic) frequencies. 1000 uniformly distributed points of Z(uS) have been used to built the Pade approximant. 
Parameters: T = 1.4 and p — 1. 


The damping [2j [33] 


7 f(k) = -\ojf(k) 


^o(fc) 

S(k) 


(8) 


For C t the following expression have been used [2 [33] : 
k 2 TG(k) 


C t 


„ (t) +(*)■ 


(9) 


where 


_i {Go/v) -2k 2 (G(k)/p-T) 2 k 2 G{k) 


V) = 


k 2 rQ +1 

1 


Go = ^pr 0 “ E + P T > 


( 10 ) 

( 11 ) 


G(k) = 


p (jw 74 (-tSS ? 1 + jSF +1)) 


k 2 


( 12 ) 


where r] is viscosity (we take it from MD simulations). 

For the dynamic structure factor we have used the ap¬ 
proximation following Refs. [[34] [35]: 


F s {k, t) 


exp 



(13) 


where c ~ 1 /uje and D is the diffusion coefficient (we 
take it from MD simulations). This approximation has 
perfect analytical form in the Fourier space: 


S s (k,(jj) = - 


Dk 2 /ujd 


7T -y/a; 2 ± (Dk 2 ) 2 

Dk 2 \ _ \ ^Ju 2 ± {Dk 2 ) 2 


exp 


idD 


K x 


idD 


(14) 


Here K\ is the bessel function. 


For p — 1 and T = 1.4 (LJ fluid) we find the best fit of 
Z(u) using the Viscoelastic approximation. The results 
are shown in Fig. [6] The best fit parameters used in Fig. [6] 
are in fact very close to those estimated independently 
from molecular dynamic simulations. The black curve in 

(a) is VAF obtained within MD simulation while the red 
curve in (a) shows VAF obtained within the viscoelas¬ 
tic model. Analytical continuation of VAF in viscoelas¬ 
tic model into complex frequencies is shown in (b), (c), 
(d) and (e). As follows from (a) the difference between 
viscoelastic and exact results for Z(t) is not very large 
however no brunch cut is seen in viscoelastic Z(t), see 

(b) and (c), contrary to exact Z(t). Similar result we see 
for for other densities and temperatures for LJ fluid. 

Analytical continuation given in Fig. [6] shows that 
Z(u) defined by a quite involved integral in the vis¬ 
coelastic model can be unexpectedly well approximated 
just by the analytical function with the poles at oj\ = 
±0.57 ± i0.97 and cj 2 = 1.74 ± il.87: 


Z M « E 

i=1,2 


Ai 


—ioo ± 


M 2 




(15) 


where Ai are the real adjusting parameters. Only at very 
small (hydrodynamic) frequencies this simple approxi¬ 
mation becomes incorrect. It should be also noted that 
the longitudinal part of current fluctuations, see Eq. 
gives the main contribution to Z(cd) except the peak at 
small frequencies where transverse fluctuations dominate. 

Summarizing this section, one can again conclude that 
there is no analytical approach to calculate VAF with 
enough accuracy to build analytical continuation in com¬ 
plex frequency plain. So the only alternative is the nu¬ 
merical methods described below. 
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FIG. 7. Analytical continuation of oscillator power spectrum 
from the real axis to the complex plain by multipoint Pade 
approximant: (a) shows the absolute value while (b) is the 
argument. 



FIG. 8. (a) We add gaussian noise with zero mean and 

a = 5 X 10 -4 to the oscillator power spectrum. Analytical 
continuation is shown in (b) and (c). The “main” poles show 
high degree of resistivity to noise. 


V. METHODS 


A. Pade approximation: Numerical multipont 
continued fraction algorithm 


B. Pade approximation: Illustrating test-examples 


1. Oscillator power spectrum amplitude: approximation of 
the analytical function with 4 poles. 


As the first test example we take the function 


Here we discuss the construction of the Pade approx- 
imants that interpolate a function given N knot points. 
Pade-approximants are the rational functions (ratio of 
two polinomials). A rational function can be represented 
by a continued fraction. Typically the continued fraction 
expansion for a given function approximates the function 
better than its series expansion. 

Algorithm: for a function f(x{) = Ui with values ui at 
N knots Xi, i = 1 , 2 ,3,..., A, the Pade approximant is 


Cn(x) = 


CLl 


a,2{x—x\) 


d 3 (x- x 2 ) _ 

° 4 (*-* 3 > - + 1 

N [x-x iV-i J + l 


+1 


+ 1 


(16) 


where a, we determine using the condition, C'/v ( x t ) = u^, 
which is fulfilled if a* satisfy the recursion relation 


ai=9i{xi), gi{xi)=Ui, i = 1,2,3,..., IV. 

(17) 

g (x) = i(* P -i) -g P -i(aO > 2 . (18) 

[x - x p -i)g p -i(x) 


fosc (^) / 2 2^ i ( \2 ’ (1^) 

0 v - u o) + 

This function is proportional to the oscillator power spec¬ 
trum. We take uj = 3 and 7 = 2 and build the Pade 
approximant using 300 uniformly distributed knots at 
uj G (—10,10). The result of the analytical continua¬ 
tion is show in Fig. [7] We worked with double precision. 
The relative error of the analytical approximation was 
less than 10“ 10 even in the pole-regions. 



FIG. 9. Analytical continuation of the logarithm: (a) shows 
the absolute value while (b) is the argument. 




























FIG. 10. Analytical continuation of the square root function: 
(a) shows the absolute value while (b) is the argument. Ar¬ 
ray of peaks and dips in (a) represent the branch cut. Graphs 
(c) and (d) show “exact” absolute value and argument of the 
function. Figure (e) is the density plot of the relative differ¬ 
ence between the exact and Pade-approximation. The coinci¬ 
dence is perfect everywhere except the white zone where the 
functions differ because the branch cuts of the approximation 
and the “exact function” have been chosen differently. 


2. Stability of the Pade-approximation 

We add gaussian noise with zero mean and cr = 5 x 10 -4 
to the oscillator power spectrum considered above, see 
Fig. [8] Analytical continuation is shown in (b) and (c). 
The “main” poles are still clearly seen. So analytical 
continuation by Pade approximation is quite resistive to 
noise if the noise correlation length is short enough. 


3. Analytical continuation of In-function by the Pade 
approximant 

Now we illustrate how behaves singular function in the 
complex plain when we do its Pade analytical continua¬ 
tion. We take 

/lnM =ln(l + w). (20) 

We build the Pade approximant using 300 uniformly dis¬ 
tributed knots at real uj G (0,50). The result of the 
analytical continuation is shown in Fig. [9] There is cut 
(—oc,—l) in the complex cj-plain. Analytical continua¬ 
tion based on the Pade approximant reproduces the cut 
by the array of poles and knots (where f\ n = 0), see 
Fig.U Away from the cut the accuracy of the analytical 
continuation is satisfactory as in the upper illustrating 
example while |cj| < 50. 



FIG. 11. Analytical continuation of tanh(cj): (a) shows 
the absolute value while (b) is the argument. We build the 
Pade approximant using 500 uniformly distributed knots at 
u G (-10,10). 


4- Analytical continuation by the Pade approximant of the 
function with the square root singularity. 

Finally we take the function with the square root sin¬ 
gularity to test the Pade-approximation: 

/sqrM = Vu> - i. (21) 

We build the Pade approximant using 300 uniformly dis¬ 
tributed knots at real uj G (—10,10). Then we analyt¬ 
ically continue the Pade polinomial (it is in fact com¬ 
plex even at real uS) to the complex u;-plane as shown 
in Figs. [lOja) and (b). Array of peaks and dips in (a) 
represent the branch cut: this is typical for pade approx¬ 
imation. Graphs (c) and (d) show “exact” absolute value 
and argument of the function. The branch cut parallel 
to the real axis is typical choice for “computer” build 
in functions (we have used Mathcad). Pade approxima¬ 
tion have chosen different direction for the branch cut, 
parallel to the imaginary axis, see (a) and (b). Fig¬ 
ure (e) is the density plot of the absolute value of the 
difference between the exact and Pade-approximation 

(I/ln - /in Pade) |/(l/in| - II )• The coincidence is per- 
feet everywhere except the white zone where the func¬ 
tions differ because the branch cuts of the Pade approx¬ 
imation and the “exact function” are different. 


C. Limits of applicability of Pade approximation 

As follows from the examples, if we approximate a 
function by the Pade polinomial at certain domain at the 
real axis then the analytical continuation is more or less 
perfect at the circle in the complex plain (around that 
domain) with the radius about the length of the domain. 

The branch cuts are represented by an array of poles. 

There is a problem with the branch cuts: we can draw 
them differently in the complex plain, only edges are 
fixed. Different choice of the branch cut curve corre¬ 
sponds different analytical continuation. But the Pade 
polynomial chooses the cut curve somehow “automati¬ 
cally”: we do not well control that. So the Pade approx- 
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imation is a useful tool if one needs to identify the posi¬ 
tion and types of the singularities of the function in the 
complex plain like poles and the the branch cut edges. 
For functions without branches analytical continuation 
in unique and the Pade approximation well produces it, 
see, e. g., Fig. [Til • 


VI. CONCLUSSIONS 

Singularities of dynamic correlation functions in com¬ 
plex plane usually correspond to the collective or local¬ 
ized modes. We have found that instead of few number 
of isolated poles velocity autocorrelation function of LJ 
particle system in complex plain shows the singularity 
manifold forming branch cuts that suggests LJ velocity 
autocorrelation function is a multiple-valued function of 


complex frequency. The brunch cuts are separated from 
the real axes by the well-defined “gap”. The brunch cuts 
are quite stable with the respect to temperature and den¬ 
sity variation. We have found the trace of the singularity 
gap at temperatures several orders of magnitude higher 
than the melting temperature. Our working hypothesis 
is that the branch cut origin is related to the interference 
of short-time one-particle kinetics and and long-time col¬ 
lective motion (hydrodynamics). 
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